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Abstract 

The monarch butterfly annually migrates from central Mexico to southern Canada. 
During recent decades, its population has been reduced due to human interaction with 
their habitat. We examine the effect of herbicide usage on the monarch butterfly's 
population by creating a system of linear and non-linear ordinary differential equations 
that describe the interaction between the monarch's population and its environment 
at various stages of migration: spring migration, summer loitering, and fall migration. 
The model has various stages that are used to describe the dynamics of the monarch 
butterfly population over multiple generations. In Stage 1, we propose a system of 
coupled ordinary differential equations that model the populations of the monarch 
butterflies and larvae during spring migration. In Stage 2, we propose a predator-prey 
model with age structure to model the population dynamics at the summer breeding 
site. In Stages 3 and 4, we propose exponential decay functions to model the monarch 
butterfly's fall migration to central Mexico and their time at the overwintering site. 
The model is used to analyze the long-term behavior of the monarch butterflies through 
numerical analysis, given data available in the research literature. 



1 Introduction 



The migration of the monarch butterfly is a marvel of nature. It is a journey through 
time and space that spans a distance of 4500 km and at least four generations of monarch 



butterfly, the exact number of generations depends on local climatological conditions 12 



The persistence of this yearly cycle is heavily dependent on two plants: the milkweed 
plants, of the family (Asclepiadaceae) found all over North America, and the Oyamel 



fir tree found at its overwintering habitat 16 . The milkweed is the only food source 
for the larvae and the Oyamel fir trees help keep the monarch in a cool state during its 
winter hibernation |12| . Due to deforestation in the mountains of central Mexico and the 
increased usage of herbicide in the United States and Canada, the plants the monarch 
butterfly depends on have been reduced [5]. 

We are primarily concerned with the effect herbicide usage has on the long-term pop- 
ulation dynamics of the monarch butterfly. For this purpose, we develop a multi-stage 
model that describes the monarch butterflies migration and use this information to create 
a discrete time model of the behavior of the population year after year. Our model is 
based on values obtained from previous research on the monarch butterfly migration and 
the milkweed plant. We then use our model to estimate the impact of herbicide usage on 
the population of monarch butterflies. 



2 Biological Background 
2.1 Monarch Butterfly Life-Cycle 

The monarch butterfly {Danaus plexippus) is rare among migratory animals and unique 
among insects. In the family of insects, the desert locust is the only other species that mi- 
grates a comparable distance |15| . Desert locusts have a dynamic migratory cycle, a cycle 
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Life Cycle of the Monarch Butterfly: 

Egg Larva 




Adult Butterflv Pupa 
Inn ages from http://www.rn onarch-butterflv.com/ 



Figure 1: The life-cycle of the monarch butterfly has four stages: egg, larva, pupa and 
adult 

dependent on "directed movement controlled by tides or wind, with navigation abilities 



not essential" 15 , unlike the monarch butterflies, which have a seasonal migration. It is 



rare among migratory animals, because the generation that leaves the overwintering site 
in central Mexico, in the spring, is not the generation that returns to the overwintering 



site the following fall 16 . There are multiple subspecies of monarchs, migratory and 



nonmigratory [14]. A subspecies known as Danaus plexippus plexippus are the migratory 
monarchs, which will be the primary focus of this paper. There are several migratory 
populations of monarch butterflies as well; migrant monarch butterflies that live east and 
west of the Rocky Mountain range. We focus primarily on the populations east of the 
Rocky Mountain range, because they have the largest population and the longest migra- 
tion route [18] . 

There are several stages of the monarch butterflies flight: the spring migration from 
central Mexico to southern Canada, the summer loitering in southern Canada, and the fall 
migration from southern Canada back to central Mexico [13] . We will follow the convention 
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of Lincoln P. Brower and designate the monarch populations that are traveling, either from 
south to north or from north to south, as migrants |3j. 

The monarch butterfly migration begins in the Oyamel fir trees on the mountain of 
Sierra Palon in central Mexico [3], where they spend the winter in a state of torpor 
and reproductive diapause, a state of non- reproduction [18]. Research has shown that 
shorter day length, lower temperature, and larvae feeding on older milkweed increases 
the likelihood that a monarch will enter reproductive diapause, the state necessary for 
fall migration ml. Though the mechanisms are known, the exact cause of the monarch 
butterfly's transition to its autumnal migratory state is unknown and is a current subject 
of research [13| . 

When the unknown mechanism is triggered, the monarch butterflies begin their migra- 
tion, leave reproductive diapause, and become reproductively active pi. While reproduc- 
tively active, the female monarch can lay up to 700 eggs during her lifespan of seven to nine 



months 16 . Female monarchs search for young milkweed leaves to lay their eggs, laying 
one egg per milkweed leaf, before flying off to find another milkweed plant to lay more 
eggs [ 16] , After laying an egg, the female monarch resumes her flight north, continually 
laying eggs until she dies. 

Meanwhile, the offspring hatch from their eggs after three to four days [16| . The 
larvae begin life by consuming portions of their egg before moving on to eat the milkweed 
plant [16]. The larvae have five stages of growth, called instars. The first four instars 
end after each larval molt and the final instar ends when the larvae become pupae. The 
complete larval stage lasts approximately two weeks, where the larvae spend the entire 
stage on one milkweed, during which the larva grows to about 2500 times its original 
size 



16 



The monarch larvae search for a dark place to begin their pupal stage. This stage lasts 
approximately ten days, during which the entire structure of the larvae breaks down to 
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be reconstituted into their adult form 16 . Urquhart noted that temperature can either 
retard or accelerate the growth rate at every stage of development in the monarch butterfly 
life cycle. This means that the ten days given for the monarch pupal stage, like the other 
values, are averages. 

At the end of the pupal stage, the larvae become adult monarchs ready to resume 
the migration begun by its parents. Unlike its parents, its life is reduced by a significant 
amount, living only two to six weeks, whereas its overwintering parents lived up to nine 
months [16] . We simulate this part of the life-cycle of the monarch butterfly in Stage 1 of 
our model. The migration continues in this fashion, parents beget larvae, the parents die, 
the larvae grow up and fly further north. 

Lincoln P. Brower determined the geographical extent of each generation of monarch 
butterflies through chromatography analysis of the cardenolides, the toxic chemical found 
in milkweed plants, inside each monarch butterfly 13] . The cardenolides in different species 
of milkweed plants have specific chemical profiles and each of these milkweed species is 
located within different geographic ranges. Brower used these two facts to determine that 
"the first spring generation is produced largely in Texas and Louisiana" and "continue the 
migration northwards to the Great Lakes region and Southern Canada" |4j. 

We model the next phase of the migration, the time the monarch stays in Southern 
Canada, in Stage 2 of our model. This spring migration usually begins in the middle of 
March and ends in early June |16| . 

The monarch butterflies continue their life-cycle in Southern Canada and the Northern 
United States. This stage usually lasts from mid-June to mid-August. At the end of 
this stage, the monarch receives environmental cues that cause it to enter reproductive 
diapause [16]. The monarch butterfly becomes a migrant and begins its return south 
toward the overwintering sites of the previous fall. We simulate this behavior in Stage 
3 of our model. When the monarch enters reproductive diapause, it increases its lipid 
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stores by constantly feeding on nectar El. The monarch butterfly needs this lipid reserve 
to survive the winter, during which it will feed at a much reduced rate |4j. Once migrants 
arrive to central Mexico they enter a state of hibernation which we simulate with Stage 4 
of our model. Unlike the spring migration which is composed of multiple generations of 
migrants, the fall migration consists of only one generation of migrants |18| . 

2.2 Common Milkweed (Asclepiadaceae syriaca) Life-Cycle 

An understanding of the monarch butterfly life-cycle would be incomplete without some 
discussion of the milkweed family (Asclepiadaceae) of plants. Milkweed plants produce 
cardenolides (cardiac glycosides) and latex [I] as their primary defensive measures. We 
primarily focus on the life-cycle and development of the common milkweed plant and the 
effect of herbicide on its development. The common milkweed is a perennial plant that 
reproduces primarily from shoots off the main plant colony, but also reproduces from 
seeds [9j. It typically grows between 1 and 1.5 meters tall [jjj. The young leaves are the 
preferred site for the monarch female to lay her eggs [16] . 

The cardiac glycosides found in the common milkweed are toxic to many animals, 
because it "can disrupt the ionic balance of a number of different cell types in animals, 
including heart muscle, vascular smooth muscle, neurons, and kidney tubules" (9l. An 
indication of the level of destruction of the milkweed by herbicide can be found in Iowa. 
In 1999, the common milkweed was present in approximately 50% of Iowa corn and soybean 
fields. In 2009, the percentage of common milkweed was present in only 8% of the fields [8j. 
Since much of the Midwest is farmland, this is an indication that a significant portion of 
the monarch butterfly habitat is at risk. 
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2.3 Previous Mathematical Work 

Since 1960, Lincoln Brower, Fred Urquhart, other zoologists, and other biologists have 
tried to understand the life-cycle and migration of the monarch butterfly. Despite all the 
work conducted on monarch butterflies from a zoological and biological perspective, there 
has been a dearth of mathematical work. To our knowledge the only other work was a 
discrete model by Yakubu et al. [17j. The authors modeled the migration and life cycle of 



the monarch butterfly with a set of difference equations 17 . They assume that the spring 
migration, from the overwintering site in Mexico to Southern Canada, consists of three 
generations and the fall migration, which travels from southern Canada back to Mexico and 
consists of one generation. This comprises a total of 4 generations throughout the entire 
cycle. The goal of their project was to investigate the persistence of the monarch butterfly 
population, with a spatially discrete advection model with emphasis on compensatory 
(contest competition) and overcompensatory (scramble competition) dynamics |17| . 

In their work, the authors assumed non-stochastic extinction of the population and 
discrete reproduction during the spring migration. A threshold parameter, or basic re- 
productive number, for the persistence or extinction of the monarch butterfly population 
was found and they analyzed it in different situations. Based on their findings, extinction 
or persistence of the population in generation 4 depends on the non-migratory popula- 
tion size in generation 3. Different behavior is observed with different parameters and 
non-migratory population size. 

In our model, we assume that during their migration, the monarch butterflies reproduce 
continuously. We also model their movement by advection. Finally, we include the major 
larval food source, the common milkweed. We numerically investigate the long-range 
behavior (over 30 years), of our model. 
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3 A Model Under Consideration 



3.1 Description of the Model 

The difficulty of modeling the monarch butterfly arises from its unique migratory nature. 



We analyzed data from the website Journey North 11 and found that monarch butterfly 



stays within a temperature range of approximately 15.5 °C to 24 °C, see Figure 2(b) 



We also analyze the first monarch sightings of the year, available on the website Journey 
North [111, and found that monarch butterflies travel at an approximately linear rate, see 



Figure 2(a) With this, we assume that the amount of milkweed available to the monarch 
butterfly is constant, because the monarch butterfly constantly moves north into areas of 
new milkweed. Thought of in another way, we assume that the monarch butterflies occupy 
an "expanding box" traveling at a constant speed north. 



y = 0.13B«*« +2045 



Migration of the Monarch Butterfly 
MP* 
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(a) 2011 monarch butterfly spring migration (b) Temperature Data for First Monarch Sightings 

by Day (2010) 

Figure 2: 



We divide the model of the migration cycle into four stages. The first stage focuses 
on the monarch migrants' flight from central Mexico to southern Canada. This stage 
incorporates multiple generations, since the monarch butterfly continually reproduce along 
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the way. We consider the initial generation as generation zero, those monarch butterflies 
that survive the hibernation phase in central Mexico. We consider generation i as the 
monarch butterfly that arrive in southern Canada (according to existing data over the 
years, 3 < i < 7). 

The second stage describes the butterfly in southern Canada and takes into consider- 
ation the effect of herbicide on the common milkweed. We assume an age-structure on 
the monarch butterfly, larva and adults. We also assume a predator-prey model, with the 
larvae as the predator and the milkweed as the prey. 

The third stage of the model simulates the monarch migrants' return to central Mexico 
from southern Canada. We model this with an exponential decay function, because the 
monarch butterfly is non-reproductive during the fall migration, yet continually die due 
to various environmental factors (weather, natural catastrophes, etc. ... [7]). 

The fourth stage models the hibernation phase and we also simulate this stage with 
an exponential decay function, because the monarch butterflies die during their dormant 
phase. 

3.2 Stage 1 

Our model for the first stage is presented in the following form and the parameters are 
shown in Table [U 

We propose the following system of equations to model the population dynamics of 
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Figure 3: Flow of monarch butterfly reproduction from first to last generation in Stage 1. 

the monarch butterfly during the spring migration: 

dM 



dt 
dLt 

~dT 
dM x 

dt 
dL 2 

dt 
dM 2 
dt 



dLj 

dt 
dMi 

~dT 

with initial conditions 



-»qM , (1) 

aiMoA)- ( 7 + /xi)Li, (2) 

jLx-^Mx, (3) 

aiMxAo- (7 + /xi)L 2 , (4) 

7L2 - M2-M2, (5) 

aiMi_iA - (7 + L h (6) 

jLi-^Mi, (7) 



Mo(to) / Ljito) = Mj(t Q ) = Ofor 1 < j < i, 
and with the parameters given in Table [T] In Equation [TJ Mq represents the monarch 
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Table 1: Parameters for Stage 1 (for more details see Appendix |A|) 
Parameter Biological Meaning Default Value Unit 



fiQ Death rate of overwintering monarchs 0.1198 - 0.1997 1/day 

a\ Growth rate of larvae 2.232 1/day 

Aq Percentage milkweeds not killed by herbicide 1 

7 Maturation rate of larvae 0.03571 1/day 

/xi Death rate of larvae 0.0902 - 0.1397 1/day 

fj,2 Death rate of adult monarchs 0.07143 1/day 



butterfly population in central Mexico. In Equations [2] and [3j L\ and M\ represent the 
population of the first generation larvae and monarch butterflies who migrate to Canada 
in the spring. In Equations [4] and [5j L<i and M2 represent the populations of larvae 
and monarch butterflies, in the second generation who migrate to Canada, respectively. 
Similarly in Equations [6] and YR from the model L« and Mj are the population of larvae 
and monarch butterflies in the i th generation. 

In Stage one Equation[TJ the term — HqMq describes the rate of change in the population 
of fall migrants, the monarch butterflies that overwinter in central Mexico. The term ^lq is 
the death rate of the initial generation of monarch butterflies. These monarch butterflies 
were in reproductive diapause and they will not produce monarch butterflies who are in 
reproductive diapause, so the population decreases at a proportional rate. 

In Equation [2j the parameter a\ is the per capita growth rate of the larvae, given in 
Table [T] The rate a\ is multiplied by the initial population leaving Mexico and the amount 
of milkweed, which is a ratio, along the way since the larvae depend on the milkweed. 
Thus, the amount of larvae, L\, is proportional to the milkweed, Aq. The second term 7 
in Equation [2] is the maturation rate of the monarch butterfly, from larva to adult, so we 
have the term — 7L1. We also consider those larvae that die before they become adults. 
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They leave the population, so we have the term —^iL\. 

Next, we look at Equation [3] in stage one. 7 is the maturation rate of larvae. The 
monarch butterflies also die at a rate proportional to their population, thus we have — fj,\ 
multiplied by Mi, the number of monarch butterflies in generation one. This pattern 
continues until the i th generation (the total number of reproduction generations in stage 
one), when the monarch butterfly population reaches southern Canada. 

3.3 Stage 2 

We model stage two with the following system of equations: 



= -(7 + Ml) L m, (8) 

= -HiM in + jL in , (9) 
d ^ = a 2 AM m -( 1 + f i 1 )L s , (10) 
^= 7 L S -^ 3 M S , (11) 

A A / A \ 

-^ = aA(l--)-A(a + (3L s ). (12) 

The previous equations have the initial conditions: 

A(h)^0, L s (h) = M s (h) = 0, M in {h) = Mi(h), L in (t 1 ) = L i {t 1 ). 

and the parameters listed in Table [2] In Stage 2, we consider the interaction between the 
larvae, the milkweed, and the adult monarch butterflies and also consider the effect of 
herbicide on the milkweed. We obtain the visual illustration of the model in Figure [4} 

The terms Lj n and Mj n are the population of larvae and monarch of the last generation 
of stage one, respectively. The terms L s and M s are the population of larvae and monarch 
adults that are in reproductive diapause. A small portion of the larvae mature and become 
adult butterflies and this is reflected by the term —^Li n in Equation [8l where 7 is the 
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Para. 


Table 2: Parameters for Stage 2 (for more 
Biological Meaning 


details see Appendix 
Parameter Value 




Unit 




7 


Maturation rate of Larvae 


0.03571 


1/day 




Ml 


Death rate of Larvae 


0.0902 - 0.1397 


1/day 




M2 


Death rate of later monarchs 


0.07143 


1/day 




Q?2 


Growth rate of larvae with Monarchs 


2.6 


m 2 /(kg 


* day) 


M3 


Death rate of non-reproductive monarchs 


0.005 


1/day 




a 


Growth rate of milkweed 


0.007 


1/day 




K 


Carrying capacity of milkweeds 


1.79188 


kg/m 2 




a 


Percentage of milkweed destroyed by herbicide 1 


1/day 




P 


Consumption of milkweed by larvae 


5 • 10~ 9 


1 / (larvae * day) 



maturation rate, see Table [2] The parameter ^ is the mortality rate of the adult monarch 
butterfly. Equation [9] represents the population increase of reproductively active monarch 
butterflies. 

The larvae of Mj ra are denoted by L s and they either die or they mature, becoming 



M s . We model this with Equations 10 and [IT] of Stage two. Equation 12, describes 
the interaction between the milkweed, the larvae, and the herbicide. This system is a 
modified predator-prey system where the larvae are the predator and the milkweed is the 
prey. Adding the adult butterfly to the system changes the dynamics of the system, since 
the butterfly has a positive effect on the plant through pollination [16]. Since there are 
multiple pollinators of the milkweed, we consider the effect of the adult butterfly negligible. 

3.4 Stage 3 

In the third stage, the butterflies are migrating back to central Mexico and are in reproduc- 
tive diapause, therefore their populations doesn't increase. Also, a significant proportion 
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Figure 4: Interaction between larvae, adult, and milkweed. 

of them do not reach central Mexico due to the following factors: natural catastrophes, 
weather, and other environmental factors [7]. Thus, we have an exponential decay function 
of the following form, where \iu n represent the death rate at this stage. The number of 
butterflies in reproductive diapause that leave southern Canada, M s {t2), are the initial 
condition of the following equation: 

dM fin 



dt 



-fifi n Mfi n . (13) 



3.5 Stage 4 



The last stage models the butterfly in its dormant phase. This is qualitatively similar to 

Stage 3, because they are still in reproductive diapause, and a large portion of them will 

die. We model this with an exponential decay function using a different rate than Stage 

3. The population of butterfly that arrive in central Mexico, from stage three, represent 

the initial condition, M w (t^), of the following equation: 

dM w , . 

= -fi w M w . (14) 
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Table 3: Parameters for Stage 3 and 4 (for more details see appendix [C]) 
Parameter Biological Meaning Default Value Unit 

Hfin Death rate of migrant monarchs 0.0056 1/day 

jj, w Death rate of overwintering monarchs 0.0042 1/day 

4 Numerical Results 
4.1 Preliminary Results 

From the data we obtain from the Journey North [II] , we see a correlation in the distri- 
bution of milkweed and monarch butterfly as they travel from the overwintering site in 
Mexico to southern Canada. 



Figure 5(a) shows the distribution of both new leaves and adult butterflies through- 
out the United States, with respect to different latitudes (from the southern to northern 
sections of United States). Again with data from Journey North [II], we obtain a similar 
illustration for other monarch butterfly sightings for both the fall and spring migrations, 
the complete annual migration, as shown in Figure [5 (b)| 

4.2 One Year 

In this section, we look at the behavior of the monarch butterfly population over one year, 
under various conditions. We start with an initial monarch population of 150, 000, 000, 
based on the approximation found on the Journey North website [11] , an initial milkweed 
percentage of 0.6 and 0.4, and we use the various parameters in Tables [I] through [3] We 
generate the annual population behavior shown in Figures [6(a)] and [6(b)] 



Next, we run the simulation with different values of a, the herbicidal rate. As we see in 
Figures 6(a) through 7(d)| we have different behavior when we vary the value of a. First, 



we consider Aq (initial milkweed percentage in kg per square meter) equal to 0.4 and a 
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(a) 2011 monarch butterfly spring migration (b) 2010 monarch butterfly fall and spring migra- 
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Figure 5: First Monarch Sightings (2010) [IT] 

equal to 10 and we obtain the graphs in Figure [6j 

When simulate our model over a year we observe unique dynamic behaviors. The 
generations of the monarch butterflies overlap grow over time. When they arrive at Canada 
we observe a different behavior. There is a decline in the monarch and larvae population. 
When this simulation is repeated over many years the butterfly population at the end of 
each cycle exhibits discrete logistic-type behavior. 

If the constants Aq and A(ti) = Ai are sufficiently small, we observe a change in 
behavior. The graphs of the population, with respect to time, are inverted. This suggests 
a critical point in the graphs. When the simulations are run close to the turning point we 
note that the range of values becomes decreases. We found the critical point numerically 
and it's value is approximately Aq = Ai = 4.019. When looking at this over one year, it 
we see that a higher value of A will show logistic growth and a lower value of A will show 
logistic decay, depending on how extensive the herbicide is used. 

When we vary a in the equations we observe that if a is small, then the population 
is more localized in one area. In fact, increasing a by a factor of 10 eventually leads to 
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(a)Annual monarch butterfly population behavior (b)Annual monarch butterfly population behavior 
in the United States with Ao = 0.6 and cr = 1 per in the United States with Ao = 0.3 and <j = 1 per 
day. day. 
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(c)Simulation of the monarch butterfly population (d)simulation of the monarch butterfly population 
in the United States over a 30 year period with in the United States over a 30 year period with 
Aq — 0.6 and a — 1 per day. Aq — 0.3 and a = 1 per day. 

Figure 6: Variation of the amount of milkweed along the migratory route. 
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(a)Annual monarch butterfly population behavior (b)Annual monarch butterfly population behavior 
in the United States with Aq = 0.4 and a — 10 per in the United States with Aq = 0.4 and <j = 1 per 
day. day. 
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(c)Simulation of the monarch butterfly population (d)Simulation of the monarch butterfly population 
in the United States over a 30 year period with in the United States over a 30 year period with 
Aq — 0.4 and a = 10 per day. Aq — 0.4 and a = 1 per day. 



Figure 7: Variation of a for various values of A\. 
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(a)Annual monarch butterfly population behavior (b)Prediction of the monarch butterfly population 
in the United States with Aq — and a = 100000 in the United States over a 30 year period with Aq = 
per day. and a = 100000 

Figure 8: Large values of a lead to extinction. 

extinction, while smaller values of a shows stabilization. 

To verify our model reflects reality, we simulated the extinction of the milkweed and in 
our model the monarch butterfly population dies off as expected. If we simulate the pop- 
ulation of milkweed without harvesting from herbicide, we see that the monarch butterfly 
population increases by a factor of 2000. 

We vary the amount of milkweed at the various stages of the model, in the United 
States and Canada, in our simulations. We see the monarch butterfly population is more 
sensitive to milkweed in the southern United States than other areas. 

For Aq = 1, and Ai = 1 there is not much difference when a = 0,1,10 or 100 but 
there is a difference when a = 1000. For instance, the first four values of a, the value 
approaches 10 10 on a logarithmic scale and when a = 1000 it approaches 0. We note that 
when a = 100 it has a moderate effect on the population. 

For Aq = 2 and Aj = 2 the graphs are similar, i.e. they have approximately the same 
upper horizontal asymptote, but when a varies the only change is the shape. 
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(a)Annual monarch butterfly population behavior (b)Prediction of the monarch butterfly population 
in the United States with Aq — and a = 100000 in the United States over a 30 year period with Aq = 
per day 1.78 and a — 

Figure 9: Lack of herbicidal spraying leads to population stability. 

For Aq = 0.2 and A, = 0.2 we see when a = the curve approaches a nonzero 
asymptote but when a = 1, the curve approaches a zero asymptote, and when a = 10, 
the population quickly dies off and so on with higher values of a. The curves are more 
sensitive to higher value of cr's. 

For Aq = 0.5 and A{ = 0.5, we observe that for a = 0, the curve has a large asymptotic 
value, but when a = 1 this asymptote is lower. When a = 10 the curve decrease to about 
third of its initial value. When a = 100 the value approaches zero. In the simulation, 
we see that for a = 1000 it has oscillatory behavior. This graph also shows another 
asymptote, which the monarch butterfly population tends to go to a = 10. 

For Aq = 0.5 and Ai = 0.4 and a = the population converges to an asymptotic 
line. For a = 1000, the population approaches extinction. Although this level of herbicide 
is excessive it is reasonable because with the amount of herbicide used, it will cause to 
population to go to zero. 
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(a) Annual monarch buttery population behavior in 
the United States with Ao = 0.4, At = 0.7, and 
(7 = 1 per day. 
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(c)Simulation of the monarch buttery population 
in the United States over a 30 year period with 
A = 0.4, Ai = 0.7, and a = 1. 

Figure 10: Milkweed population along the i 
term population size. 
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(b) Annual monarch buttery population behavior in 
the United States with Ao = 0.7, A t = 0.4, a = 1 
per day. 
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(d) Simulation of the monarch buttery population 
in the United States over a 30 year period with 
Ao = 0.4, Ai = 0.7, and a = 1. 

tigratory route has significant effect on long 
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Figure 11: Simulation of the monarch buttery population in the United States over a 30 
year period with Aq = 1 and A{ = 1. 
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Figure 12: Simulation of the monarch buttery population in the United States over a 30 
year period with Aq = 2 and Ai = 2. 
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Figure 13: Simulation of the monarch buttery population in the United States over a 30 
year period with Aq = 0.5 and Aj = 0.5. 
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Figure 14: Simulation of the monarch buttery population in the United States over a 30 
year period with Aq = 0.2 and Aj = 0.2. 
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Figure 15: Simulation of the monarch buttery population in the United States over a 30 
year period with Aq = 0.5 and Aj = 0.4. 

4.3 Estimation of Population Density 

Obtaining an accurate population count of all monarch butterflies is very difficult, if not 
impossible. We attempt to find the population density at specific times, so that future 
experiments can either verify or refute our model. 

We chose three coordinate points on the map of the United States to create three 
vertices of a triangle. This triangle covers the area where monarch butterfly activity 
was observed. The three points are: A = (23N, lOOVF) at the central Mexico, B = 
(48iV, imW) at the Northwestern U.S. and C = (48iV, 71W) at the Northeastern U.S. 

In Figure [i~6[ the left side of the triangle lies along the eastern side of the Rocky 
Mountain range. We assume the area is a right triangle. We consult Google Earth, an 
online geographical information system, to find the area of the triangle, approximately 
3,111,352 km 2 . 
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Figure 16: An approximation of the area or reported monarch butterfly activity. 

In order to calculate the density for each day, we divide the area into 75 horizontal 
segments with equal height, where each segment represents the distance traveled by the 
monarch butterflies in one day of the spring migration. This is not a robust assumption, 
because the population is likely distributed in multiple strips. 

We index the numbers from 1 through 75, from South to North. Based on the formula 
for the area of a right triangle, the formula for calculating the area of a particular strip is 
given as: 

Area^ = • (* 2 - (i - I) 2 ) 

We presume the entire population of monarch butterflies in the i th strip at the i th day 
of the spring migration. The 75 th strip represents the second stage. The (75 — i) th strip 
is at the i th day of the fall migration and finally the I s * strip represents the forth stage. 



We calculate the population density at each specific time as shown in Figure 17 



5 Discussion 

If we simulate the model over the course of a year, we observe unique dynamic behavior. 
The monarch butterfly generations overlap one another and the generations grow with 
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Days From March 15 to the Next Yea 



Figure 17: population density of the monarch butterfly over four stages of migration 

respect to time. When the monarch butterflies arrive at Canada, we observe different 
behavior, there is a decline in the monarch and larvae population. This is the behavior 
we observe in nature, because the monarch butterfly is stationary in Canada and they 
stay for an extended period of time (limited food resorces). This means the milkweed 
population changes with respect to time. When we repeat the simulation over many years 
the monarch butterfly population shows logistic growth. 

If we set the constants Aq and ^4(100) = A{ at a low enough value, then we see a change 
in behavior. The graphs of the population invert with respect to time. This suggests a 
critical point. When we run simulations closer and closer to the critical point, we notice 
the range of values decrease. When we look at this for one year, we see that higher values 
of A exhibit logistic growth and lower values of A exhibit logistic decay depending on the 
value of a. 

A sensitivity analysis of a shows that for small a the population is localized in one 
area. Increasing a by a factor of 10 will lead to extinction while smaller values of a show 
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stabilization. 

We model the population when there is total extinction of the milkweed, to verify our 
model behaves as expected and the monarch butterfly population becomes extinct as well. 
We also look at the effect a value of a = has on the butterfly. We see that the population 
increases by a factor of 2000, as time increases. 

We vary the values of the milkweed in the United States and in Canada through our 
simulations. We see a larger impact on the overall population if there is more milkweed 
present in the United States than in Canada. This difference is maybe due to a larger 
supply of milkweed during their migration and the larger portion of time spent in the 
United States. 

This value changes as we change the value of a. We conclude that herbicide has a large 
effect and a reduction of herbicidal spraying is needed to stabilize the monarch butterfly 
population. In 2002, a severe winter storm in central Mexico caused the death of approx- 
imately 80% of the monarch butterfly population, at the central Mexico overwintering 
site [11] . The population rebounded the next year, though [Tl] . There was likely enough 
milkweed for the monarch butterfly population to increase. In later years it seems that 
the monarch butterfly population is smaller than the average value over the past several 
years |11| , indicating it may converge to a state. Oscillating monarch population is due 
much to detrimental weather and declining forest population and we can only say that the 
herbicide has an effect on where this oscillation should be. 

6 Further Work 

For a more accurate picture of these simulations, more data collection sites located south of 
Cape May, NJ and west of Chincoteague, VA would be instructional. These collection sites 
would be valuable because we know that during the fall migration, the monarch butterfly 
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travels in a southwestern direction. This means the data at all of these sites could be 
used to find an approximation of the population via method of numerical integration. It 
would also be advantageous if Chincoteague revived their data collection simultaneous to 
the other site. Not only should these three sites be operating at the same time, they 
should all be using the same method of collection. Further research needs to be done on 
the population dynamics of milkweed to obtain a more accurate model of the butterfly 
population. Another potential research topic related to dynamics is the effect of the 
Oyamel fir forest, in Michoacan, Mexico, on the monarch butterfly population. 
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Appendix A: Parameters for Stage 1 

The monarch butterfly start out in central Mexico at a latitude of 20° on the 75 th day of 
the year. For calculation purposes we re-scale time by setting 75 to day 0. Because the 
eggs are laid in the southern United States (30° — 35°) we can assume that 5% of the adult 
monarchs will remain in this area. Referring to the above figure the butterflies will reach 
35° approximately 25 days after departure. The general solution to the first equation in 
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our model is: 

M = M 0l ■ e~ wt 
Then to estimate our parameters we obtain: 

0.05A% = M 0i e- 2 ^° 

by substitution 

^ = -log(0.05)/25 pa 0.1198 



According to the Figure [5(a) the butterflies will reach 30° latitude 15 days after departure. 
Then using the same approach as above we obtain jUo = 0.1997. This gives us a range of 
values for 

According to Urquhart, the average time for a larva to mature from egg to butterfly 



is approximately 28 days 16 . This means there is a maturation rate 7 = 1/28 per day. 
The parameter \i\ is calculated from the high mortality rates of the larvae. From Dively 
et al 92 to 98 percent of the larvae do not make it to the adult stage. Since it takes 28 
days for the larva to mature 92 to 98 percent of the larvae population will die within 28 
days. Then using the general exponential solution we obtain values from 0.0902 to 0.1397 
per time for our 

A monarch female lays from 500 to 700 eggs over her lifespan [16] and her average 
life-span is approximately four weeks, or 28 days [16] . This means that each each female 
monarch lays, on average, between 17 and 25 eggs each day, approximately 17.857 to 25 
eggs. This gives a reproductive rate of 8.929 to 12.5 per monarch butterfly. Not all eggs 
are entirely fertile and the number of infertile eggs can be as high as 55% [10] . This gives 
a range of 4.91 and 12.5 larvae per monarch, which we use for our values of a\. 

In stage 1, the adult monarchs live anywhere from 2 to 6 weeks. Thus to average \ii 
we can assume the worst for these monarchs and assume that on average they live 2 weeks 
this means [12 = 1/14 per day. 
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Appendix B: Parameters for Stage 2 

For the parameters in stage 2, we use some of the same parameters from stage 1. Because 
stage 2 is a stopping point, we vary the milkweed in the Equation [TU] which makes our 
«2 different from a\. The parameter «2 is also different because the milkweed units have 

2 

changed. We assume a value of 2.6 k ^ ay for 02, because the milkweed population will be 

2 

decreasing not as a percentage but as in biomass thus c*2 = kg "^ ay . 

In stage 2 the butterflies in reproductive diapause are no longer mating and they 
are preparing for the upcoming journey by consuming more food. The butterflies in the 
reproductive diapause state have not been traveling either. This means the mortality rate 
/X3 will be a small value. Thus we can assume that ^ = 0.005 per day, since I///3 = 200 
days, which is approximately 6.7 months. 

Due to insufficient data the carrying capacity of the common milkweed could not be 
found, but there did exist sufficient information on the butterfly milkweed []. According to 
the grower's guide the optimal density of milkweeds was grown in the field. This included 
a density of 43,560 plants per acre. The grower's guide also included the dry weight herb 
of 104.7 pfjj: and a dry weight root of 61.9 p ^ nt - Then the weight of the entire plant is: 

104.7 + 61.9 = 166.6 



plant 

To find carrying capacity we multiply: 

166.6 • 43560 ^ • . 1 J*_ = 1.79188 % 

plant A 4050m 2 lOOOg m 2 

The growth rate, a, of the milkweed was also obtained from the data. It appears that the 
total mass of a year one plant using similar calculations as before is 10.6 g/plant and that 
of a year two plant is 132.41g/plant. To calculate the growth rate, we use exponential 
growth once more. Let the initial condition be 10.6, then we can substitute in the new 
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values in: 

A = Aie at 
132.41 = 10.6 • e 365a by substitution 
a = 0.007 per day. 

For p we know that the value of the larvae population will be on the order of 10 8 . Because 
we know that there should not be a high decline rate of the milkweed we can assume 
(3 = 5- 10-9. 

Appendix C: Parameters for Stage 3 and 4 

Upon the return trip to central Mexico an estimation of the parameter came from the 
tagged data from the Monarch Watch website data base [IT]. The data was filtered out 
through the process of having non-dated taggings removed and considering only monarchs 
that were tagged after August 8th for any year. This is important because we need to only 
consider those in migration. Then out of the ones that we are considering, the amount 
that made it to Mexico. It was found that 84 percent of the monarchs made it. Then 
we can say 16 percent of the migrants die on their journey south. From August 15 to to 
November 1 we have a time span of about 75 days. This gives a value of fifm = 0.002325 
per day. 

Then from the Journey North website it is stated that during their stay at the over- 
wintering site approximately 15 percent of the population dies off due to predation. Their 
stay at the overwintering sight from November 1 to March 15 is about 135 days. Then by 
the same method a before we obtain as parameter estimate of fi w = 0.001204. 
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Appendix D: MATLAB CODE: testa.m 

function [x,pop] = testa(P ) N,M_0,A_i,P_2,P_3,P_4,flag) 

% The function outputs the ending population of normal adult of stage 1 and 
% of the reproductively diapaused adults of stage 2,3 and 4 as a row vector x. 
% The input includes the arguments P, P_2, p_3, P_4 for the parameters of 
% differential equations of stages 1, 2, 3 and 4 respectively, flag determines 
% weather the function outputs graph or not, its default value is 1. 

% example possible input: testa( [0 . 1198,2 . 6, .6, 1/28,0. 10147, 1/14] ,4, 150000000, ... 
% . 6 , [2 . 6 , . 0025 , . 007 , 1 . 78788552 , 1 , . 000000007] , . 002298 , . 001024 , 1) ; 

/oO i AvjEj 1/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/oA 

tspan=[0,75] ; %time span of stage 1 
y0(l)= M_0; 

for i = 2:1:2*N+1 °/ set initial condition for nth generation of larvae and adult 
y0(i)=0; 

end 

[t, y]=ode45 (©monarch, tspan, yO, [] ,P) ; °/ solve system of equations of stage 1 

function yprime=monarch(t ,y ,p) %define equations of stage 1 
mu_0=p(l) ; 
alpha_l=p(2) ; 
A_0=p(3) ; 
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gamma=p(4) ; 
mu_l=p(5) ; 
mu_2=p(6) ; 



yprime(l)=-mu_0*y(l) ; °/ define equations cyclicly 
for j=l:N 

yprime (2* j ) =alpha_l*y (2* j -1) *A_0- (mu_l+gamma) *y (2* j ) ; 
yprime (2*j+l)=gamma*y (2*j ) -mu_2*y (2*j+l) ; 

end 

yprime=transpose (yprime) ; % since yprime is defined as row vector, 

% it has to be transposed to be used as column vector 



end 



•/.knitting stage i & stage 2 7o%y„y„yoyo%%7oy.y.yoyo%%%%y.yoyo%%%%y.y.yoyo%7o7oy.y.yoyo%%%%y.yoyo7o7o7o7o /. 

for i = 1:N 7 creating a matrix that stores i columns of generations' larvae population 
HTotal(i) = y(end,2*i); 

end 

for i = 1:N 7 creating a matrix that stores i columns of generations' adult population 
mlTotal(i) = y (end, 2*i+l) ; 

end 



/QT , A(^T7 O °/ / °/ / / / / a / a /°i 



■^0/ 0/ 0/ 0/0/ 0/ 0/ 0/ 0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/ 

/o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o /o 
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zO = [sum(llTotal) ; sum(mlTotal) ;0;0;A_i] ; "/sum all the entries of each column of the matric 
tspan_2 = [0,75] ; °/„time span of stage 2 
[t_2,z]=odel5s(@monarch_2,tspan_2,z0, [] , [P,P_2] ) ; 

function zprime = monarch_2(t_2,z,p) 

mu_0=p(l) ; 
alpha_l=p(2) ; 
A_0=p(3) ; 
gamma=p(4) ; 
mu_l=p(5) ; 
mu_2=p(6) ; 

alpha_2=p(7) ; 
mu_3=p(8) ; 
a=p(9) ; 
K=p(10) ; 
sigma=p(ll) ; 
beta=p(12) ; 

zprime(l)= -(gamma+mu_l)*z(l) ; °/L_in 
zprime(2)= -mu_2*z(2)+gamma*z(l) ; %M_in 

zprime(3)=alpha_2*z(2)*z(5)-(gamma+mu_l)*z(3) ; °/ L_reproductive Diapause 
zprime(4)=gamma*z(3)-mu_3*z(4) ; °/ M_Reproductive Diapause 
zprime(5)=a*z(5)*(l-z(5)/K)-z(5)*(sigma+beta*z(3)) ; %A 
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zprime = transpose (zprime) ; 

end 

'/.knitting stage 2 & stage 3 mmmmmmmmmmmmmxxmm 

m2Total = z(end,4); %take the final value of adult population 

°/QT A PT7 oVVV / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / 
/oO 1 AvjEj o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/^ 

vO = m2Total; 
tspan_3=[0,75] ; 

[t_3,v]=ode45(@monarch_3,tspan_3,v0, [] ,P_3) ; 

function vprime = monarch_3(t_3, v,p) 
mu_f in = p(l) ; 
vprime(l) = -mu_f in*v(l) ; 

end 

% plot(t_3,v); 
% v(end) 

'/.KNITTING STAGE 3 & STAGE 4 miWICIW^^^ 
m3Total = v(end) ; 

°/QT APT? A V V V V °/ °/ °/V V V V °/ °/°/ / y / / / / / / / / / / / / / / , /V D /V'/ / , / , /V D /V'/ D /V , /V D /V*/ D / , /VV D /V'/VV /VVVVVV /V 

/oO 1 fUjEj 4/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0/0 

wO = m3Total; 
tspan_4 =[0,135] ; 

[t_4,w]=ode45(@monarch_4,tspan_4,w0, [] ,P_4) ; 
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function wprime = monarch_4(t_4,w,p) 
mu_w = p(l) ; 
wprime (1) = -mu_w*w(l) ; 

end 

7, plot(t_4,w); 

°/PT flTTTMP o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/oy 0/0/ o/o/ yyyyyyyyyyyyyo/o/o/yyyyo/o/yyyyyo/o/o/o/yyy 
/or LiU 1 1 /o/o /o/o /o/o /o/o /o/o /o/o /o/o /o/o /o/o /o/o /o/o /o/o /o/o A 

ta = tspan(2); tb = tspan_2(2); tc = tspan_3(2); td = tspan_4(2) ; 
ToStore the time period (day) of each stage 

if nargin < 8 7 setting default value of flag 
flag = 1; 

end 

for i=l:N+l 

b_total(: ,i)=y(: ,2*i-l) ; 
7oCreating a matrix containing N generations of adult population 
end 

for i=l:N 

l_total(: ,i)=y(: ,2*i) ; 
7oCreating a matrix with N generations of larvae population 
end 
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if (flag) '/making this function output graphs only if flag != 



"/first graph with all the curves of 4 stages in normal scale 
figure ; 
hold on; 

"/Performing for-loops in the if-statement again, otherwise 

"/this if-statement comes up with errors . And it will comes up with 

°/b_total undefined error if thess for-loops are only instead the 

"/if -statement . Dont know why, I guess that the scopes of variables may 

"/span differently than that in Java. 

for i=l:N+l 

curvel_l=plot(t,y(: ,2*i-l) , 'b') ; 
"/plotting each generation of larvae population 
end 

for i=l:N 

curvel_2=plot(t,y(: ,2*i) , 'g') ; 
"/plotting each generation of adult population 
end 

for i=l:N+l 

b_total(: ,i)=y(: ,2*i-l) ; 
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°/ creating a matrix containing N generations of adult population 
end 

for i=l:N 

l_total(: ,i)=y(: ,2*i) ; 
/(creating a matrix with N generations of larvae population 
end 

curvel_3=plot (t , sum(l_total ' ) , 'black' ) ; 
°/ transposing the matrix, and summing each entry in each column, making 
% it a row vector and plotting it 

plot(t_2+ta,z( : ,1) , 'black') ; 
"/.plotting population of larvae_in at stage 2 

curvel_4=plot (t , sum(b_total ' ) , 'r') ; 
"/plotting total population of monarch_in at stage 1 

plot(t_2+ta,z(: ,2) , 'r') ; 
"/plotting population of monarch_in at stage 2 

curvel_5=plot(t_2+ta,z( : ,3) , 'magenta' ) ; 
%plotting stage 2 super larvae population 

curve l_6=plot (t_2+ta,z( : ,4) , ' cyan' ) ; 
"/plotting stage 2 super monarch population 

plot (t_3+ta+tb , v , ' cyan ' ) ; 
"/plotting stage 3 super monarch population 
plot(t_4+ta+tb+tc,w, 'cyan') ; 
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/(plotting stage 4 super monarch population 



xtotal = ta+tb+tc+td; 
plot ( [0 , xtotal] , [M_0 , M_0] ,' — '); 
/(plotting a horizontal line with vertical value M_0 

plot([ta ta] , [0 max(sum(l_total'))] , '--') ; 
%plotting a vertical bar that separates stage 1&2 

plot([ta+tb ta+tb] , [0 max(sum(l_total ' ) )] , ' — ' ) ; 
%plotting a vertical bar that seperates stage 2&3 

plot( [ta+tb+tc ta+tb+tc] , [0 max(sum(l_total ' ) )] , ' — '); 
%plotting a vertical bar that seperates stage 3&4 

legend( [curvel_l,curvel_2,curvel_3,curvel_4,curvel_5,curvel_6] , 'Population of ... 
adult of each generation' , 'Population of larvae of each generation' , 'Total ... 
population of adult ', 'Total population of larvae ', 'Population of next .... 
generation larvae', 'Population of reproductive diapause generation'); 

title_l=num2str (A_i) ; 
title_2=num2str(P_2(5)) ; 
title_3=num2str(P(3)) ; 

dl=title(strcat('A_i = ' ,title_l, ' , sigma =' ,title_2, ' , A_0= ',title_3)); 
d2=xlabel( 'Years' ) ; 
d3=ylabel( 'Population' ) ; 
set(dl, 'FontSize' ,12) ; 
set (d2 , ' FontSize ' , 12) ; 
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set (d3 , ' FontSize ' , 12) ; 
figure ; 

%graph only the population of adult butterfly over 4 stages, 
hold on; 

curve2_l=plot (t , sum(b_total ' ) , 'r') ; 

plot(t_2+ta,z(: ,2) , 'r') ; 

curve2_2=plot(t_2+ta,z(: ,4) , 'black') ; 

plot (t_3+ta+tb , v , ' black ' ) ; 

plot(t_4+ta+tb+tc,w, 'black') ; 

plot([ta ta] , [0 max(sum(b_total' ))],' — ') ; 

plot([ta+tb ta+tb] , [0 max (sum (b_total '))],' — ') ; 

plot( [ta+tb+tc ta+tb+tc] , [0 max(sum(b_total '))],' — ') ; 

plot( [0,ta+tb+tc+td] , [M_0,M_0] ,' — '); 

legend ( [curve2_l , curve2_2] , 'Total Population of Adult Monarch',... 
'Population Reproductive Diapause Generation'); 

al=title(strcat('A_i = ' ,title_l, ' , sigma =' ,title_2, ' , A_0= ',title_3)); 

a2=xlabel( 'Days From March 15 to the Next Year'); 

a3=ylabel( 'Population' ) ; 

set(al, 'FontSize' ,12) ; 

set(a2, 'FontSize' ,12) ; 

set (a3 , ' FontSize ' , 12) ; 
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end 

'/.setting output of the FUNCTiDiremramramraramrafflrafflmramn 

x(l) = y (end, 2*N+1) ; pending value of normal adult of stage 1 

x(2) = z(end,4); '/.ending value of super adult of stage 2 

x(3) = v(end) ; '/.ending value of super adult of stage 3 

x(4) = w(end) ; pending value of super adult of stage 4 



'/.LATER 



t"/ v v v v v v v v v y v v v v v y v v v v v v v v v v v v v v v v v v v v v v y v v v v v v v v v v v v v v v v 

» /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o^ 



temp_l = vertcat(t,ta+t_2,t_3+ta+tb,t_4+ta+tb+tc) ; 
'/.concatenating time output 

temp_2 = vert cat ( (sum(b_total ' ) ) ' , (sum (vert cat (z( : ,4) ' , z( : , 2) ' ) ) ) ' , v, w) ; 
'/.concatenating population output 
pop = horzcat(temp_l,temp_2) ; 

'/.outputs two columns, with each contains time vector and adultpopulation vector 



end 



Appendix E: MATLAB CODE: visual.m 

a=xlsread ( ' milkweedData . xlsx ' ) ; 
'/.reads the excel file 
b=xlsread( , FirstMonarch2010 . xlsx' ) ; 
c=xlsread( ' 0therMonarch2010 . xlsx' ) ; 
d=xlsread( , stuffFor2011.xlsx') ; 
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nnnny.piots the figures in the papernrar/.n 

Xplot((a(:,4)-l)*30+a(:,3),a(:,l),'r.',(d(:,4)-l)*30+d(:,3),d(:,l),'b.'); 
al=xlabel('Day of the Year 2011'); 
a2=ylabel( 'Latitude') ; 

legend( 'First Milkweed Leaves ', 'First Adult Monarch Sighting'); 
set(al, 'FontSize' ,12) ; 
set(a2, 'FontSize' ,12) ; 

nnmr/.plots the other figure in the paper^m /. /. /. 

%plot((b(:,4)-l)*30+b(:,3),b(:,l),'b.',(c(:,4)-l)*30+c(:,3),c(:,l),'r.'); 
%al=xlabel('Day of the Year 2010'); 
%a2=ylabel(' Latitude') ; 

%legend( 'First Monarchs ',' Other Monarchs'); 
%set(al, 'FontSize' ,12) ; 
%set(a2, 'FontSize' ,12) ; 

Appendix F: MATLAB CODE: density.m 

function y=density(P,f lag) 

%This function outputs a graph of population density against time (day) in 
%period of a year. P is a row vector contains all the parameters used in 
%system of equations in the 4 stages migration. 
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A_t = 3111352; 

°/ the area of a right triangle that covers the area where activities of 
y o monarch were observed. area is in unit square km. 



%creating an array of two columns with the first column storing the date 
%and second storing the area of strip of the triangle at the respective 
%date . 

for i = 1:360; 
A(i,l)= i; 

end 

%day 1 , the area is the one that at central mexico . 
A(l,2)=A_t/75~2; 

%day 2-75, the areas are calculated by the folmula for the area of right 
%triangle . 
for i=2:75 

A(i,2)= A_t*((i/75)-2-((i-l)/75)~2); 

end 

°/ stage 2, the range of monarch activity stays at the same area 
for i = 76:150 

A(i,2)=A_t*(l-(74/75)-2) ; 

end 
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°/ stage 3, area decreases 
for i = 151:225 

A(i,2) = A((226-i) ,2) ; 

end 



°/ stage 4, staying at central mexico 
for i = 226:360 

A(i,2) = A(l,2); 

end 



P_l = P(l:6); 
N = P(7) ; 
M_0 = P(8) 
A_i = P(9) 
P_2 = P(10:15) ; 
P_3 = P(16); 
P_4 = P(17); 



[a , b] =testa (P_l , N , M_0 , A_i , P_2 , P_3 , P_4 , 0) ; 



%figure out the number of data points of a column of b; 
N = numel(b(: ,1)) ; 



for i=l:N 

%round the time to the nearest integer 
temp = ceil(b(i, 1) ) ; 
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if (temp == 0) 

area(i)=A(l,2) ; 

else 

area(i)=A(temp,2) ; 

end 

%calculating density 
D(i) = b(i,2)/area(i) ; 

end 

if (flag==l) 
figure ; 
hold on; 

plot(b(: ,1) ,D, 'black') ; 
plot([0,360] , [D(l) ,D(1)] , ' — ') ; 

al=xlabel( 'Days From March 15 to the Next Year'); 

a2=ylabel( 'Population Density (l/knT2)'); 

a3=title(' Projected Population Density'); 

plot ([75 75] , [0 max(D)] , ' — ') ; 

plot ([150 150] , [0 max(D)] ,' — '); 

plot ([225 225] , [0 max(D)] , ' — ' ) ; 

set(al, 'FontSize' ,12) ; 

set(a2, 'FontSize' ,12) ; 

set (a3 , ' FontSize ' , 12) ; 

end 
end 
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Appendix G: MATLAB CODE: endingV.m 

function [y]=endingV (Y,P,flag) 

%Y time period in unit of year for the function to simulate 
%p parameters for function testa. m 

%flag optional variable, flag = means no output graph 

"/.example possible code: endingV(30, [ [0 . 1198, 2 . 6 , . 6 , 1/28,0 . 10147, 1/14] ,4, 150000000, . 6 , [2 . 6, 

P_l = P(l:6); 
N = P(7) ; 
M_0(1) = P(8) ; 
A_i = P(9); 
P_2 = P(10:15) ; 
P_3 = P(16); 
P_4 = P(17); 

for i = 1:Y 

[x(i,:),pop] = testa(P_l,N,M_0(i) ,A_i,P_2,P_3,P_4,0) ; °/„storing the results of each run 
% of testa. m function into a row, and creating a matrix that piles these rows together 
M_0(i+1) = x(i, 4); %Storing the final adult population to be used as initial M_0 for 
7 the next run. 

end 

if (nargin<3) 



48 



flag = 1; 

end 

title_l=num2str (A_i) ; 
title_2=num2str(P_2(5)) ; 
title_3=num2str(P(3)) ; 

if (flag) 
figure ; 

"/Plotting On a logarithmic Scale 
subplot (2, 1,1) ; 

plot(x( : ,4) , '-o' ) ; "/plotting final adult population in normal scale 

a3=title (strcat ( 'Monarch populations A_i = ',title_l,', sigma = ',title_2,', ... 

A_0= ' ,title_3)); 
al=xlabel ( ' years ' ) ; 
a2=ylabel( 'population' ) ; 

set(al, 'FontSize' ,13) ; 
set(a2, 'FontSize' ,13) ; 
set(a3, 'FontSize' ,13) ; 

subplot (2, 1,2) ; 

semilogy (x( : ,4) , ' -o ' ) ; "/plotting the same thing in log scale 
bl=title ( 'Monarch Butterfly Population'); 
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b2=xlabel ( ' years ' ) ; 

b3=ylabel( 'population (logarithmic) ') ; 
set(bl, 'FontSize' ,12) ; 
set(b2, 'FontSize' ,12) ; 
set (b3 , ' FontSize ' , 12) ; 

end 

y=x(:,4); 
end 

Appendix H: MATLAB CODE: finaltemp.m 

a=xlsread('FirstMonarch2010withtemp.xlsx') ;%reads the excel file 
plot((a(: ,4)-l)*30+a(: ,3) ,a(: ,7) , 'b. ') ;°/„temp vs days 
al=xlabel('Days of the Year' ); "/Labels 
a2=ylabel (' Temperature in Celsius'); 

a3=title( 'Temperatures at First Sightings of Monarchs in 2010') 
set (al , ' FontSize ' , 12) ; 
set(a2, 'FontSize' ,12) ; 
set (a3 , ' FontSize ' , 12) ; 

Appendix I: MATLAB CODE: refine.m 

function l=ref ine (a_0 , a_i , sigma_i) 

%This function plots the long term behavior of the monarch butterfly population but 
%can plot more than one according to the parameter being considered and can 
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°/ label the parameter accordingly. 

%0ne example is ref ine( [1 , 1 , 1] , [2,2,2] , [0,0,0] 

clear figure 

[q,r] =size (a_0) ; 

Y=30;%Def ault value for the number of years 
flag=l; 

for j=l:r 

%Default set for our parameters 

P= [[0. 1198,2.6, a_0(j) ,1/28,0.10147,1/14] ,4, 150000000 , a_i (j ) , [2.6,0.0025,0.007, . . . 

1 . 78788552 , sigma_i ( j ) , . 000000007] , . 002298 , . 001024] ; 
m(j,:)=endin g V(Y,P,0); 
7„semilogy(y(: , j)) ; 

end 

hold all 
for h=l:r 

semilogy (m(h, :),'.'); °/„plots it on a logarithmic scale 
s{h}=strcat ( ' sigma=' ,num2str (sigma_i (h) ) ) ; ^Display the sigma 
%s{h}=strcat('A_0=' ,num2str (a_0(h) ) ) ;%Display milkweed in US 
%s{h}=strcat('A_i=' ,num2str(a_i(h))) ;%Display initial milkweed in 
%Canada 
legend(s) ; 

end 

cl=title( 'Long Term Trends on the Monarch Butterfly Population') 
c2=xlabel(' Years') ; 
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c3=ylabel ( 'Population' ) 
set(cl, 'FontSize' ,12) ; 
set(c2, 'FontSize' ,12) ; 
set(c3, 'FontSize' ,12) ; 



